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ABSTRACT 

Stairs, Lyne & Shemar have found that arrival time residuals from PSR B1828-11 
vary periodically with a period ~ 500 days. This behavior can be accounted for by 
precession of the radiopulsar, an interpretation that is reinforced by the detection of 
variations in its pulse profile on the same timescale. Here, we model the period residu- 
als from PSR B1828-11 in terms of precession of a triaxial rigid body. We include two 
contributions to the residuals: (i) the geometric effect, which arises because the times 
at which the pulsar emission beam points toward the observer varies with precession 
phase; (ii) the spindown contribution, which arises from any dependence of the spin- 
down torque acting on the pulsar on the angle between its spin (i~2) and magnetic (b) 
axes. We use the data to probe numerous properties of the pulsar, most notably its 
shape, and the dependence of its spindown torque on fi-b, for which we assume a sum 
of a spin-aligned component (with a weight 1 — a) and a dipolar component perpen- 
dicular to the magnetic beam axis (weight a), rather than the vacuum dipole torque 
(a = 1). We find that a variety of shapes are consistent with the residuals, with a 
slight statistical preference for a prolate star. Moreover, a range of torque possibilities 
fit the data equally well, with no strong preference for the vacuum model. In the case 
of a prolate star we find evidence for an angle-dependent spindown torque. Our re- 
sults show that the combination of geometrical and spin-down effects associated with 
precession can account for the principal features of PSR B1828-ll's timing behavior, 
without fine tuning of the parameters. 

Key words: stars: rotation - pulsars: individual: PSR B1828-11 - methods: data 
analysis 



1 INTRODUCTION 

Pulse arrival times of neutron stars can be found very accu- 
rately, which allows for the determination of the spin period 
and period derivative to very high precision. Normally, the 
time of arrival residuals which are calculated by subtract- 
ing the period and the period derivative (and in some cases 
the period second derivative) are mostly white noise. How- 
ever, residuals from a small number of rotating neutron stars 
are found to exhibit long term cyclical, but non-oscillatory, 
variations with characteristic timescales of order months 
to years (Cordes 1993). The variability may be temporary 
(e.g. the Vela pulsar during its Christmas glitch [McCulloch 
et al. 1990]) or persistent (e.g. the accreting neutron star 
Her X-l [Tannanbaum et al. 1972], the Crab pulsar [Lyne, 
Pritchard and Smith 1988], and the pulsars PSR 1642-03 



[Blaskiewicz 1992], PSR B0959-54 [D'Alessandro and Mc- 
Culloch 1997] and PSR B1828-11 [Stairs, Lyne and Shemar 
2000]). The long timescales that characterize the observed 
variations would arise naturally from precession 1 , when the 
principal axes of a body (defined through the moments of 
inertia, which we will take as Ji ^ 1% ^ ^3) revolve periodi- 
cally around the angular momentum, as viewed in an inertial 
frame. An ellipticity e = (Is — Ii)/Ii <C 1 would be expected 
to produce variations in the timing residuals of an axisym- 
metric body with a period P p = P*/e ~ 3.2 P* (sec) (10 s e) -1 
years, where P* is the rotation period. The arrival time vari- 
ations characteristic of precession would be strictly periodic, 
but not sinusoidal for a triaxial rotator. 

There are two physical causes for time of arrival resid- 
uals (At) in a precessing neutron star (Cordes 1993). One 
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1 Throughout this paper, we call this phenomenon precession, 
as has become common in the literature, although purists might 
prefer the term nutation. 
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is directly geometrical: as the rotating star precesses, the 
symmetry axis of its radiation beam crosses the plane de- 
fined by the angular momentum of the star and the direc- 
tion to the observer at times that vary periodically over the 
precession cycle. The magnitude of the variability is set by 
the amplitude of the precession, which is roughly the wob- 
ble angle, 8 (defined as the angle between the angular mo- 
mentum and the principal axis corresponding to the largest 
moment of inertia, Fig. 0. Typically 6 <C 1, and the am- 
plitude of the time of arrival residuals is Ai gco ~ #P*. In 
addition, the dependence of the spindown torque acting on 
the pulsar on the angle between its spin and magnetic axes 
produces a timing residual that can be comparable to and 
even exceed At gco . If we assume that the spindown torque 
is proportional to fi — a(Cl ■ b)b, where ft — fifi is the 
angular velocity, b is the magnetic axis, and the dimension- 
less parameter a is a measure of the angular dependence 
(a = 1 for a spinning magnetic dipole radiating into vac- 
uum) then the spindown rate varies over the precession cy- 
cle as well, producing a timing residual At s d ~ a8Pp/t s ^ ~ 
(aPp/P*t s d)At gco , where t sc j is the spindown timescale for 
the pulsar; the dimensionless parameter r s( j = Pp/P*i s d ~ 
3.2 Pp (years)/P*(sec)(i s d/10 7 years) may be large. Associ- 
ated with these arrival time residuals are period residuals 
(AP/P*) gco ~ 6P*/P P w 3.2 x 10" 8 6»P*(sec)/P p (years) and 
(AP/P*) s d ~ ar sd (AP/P*) g co. 

The best candidate to date for truly periodic long term 
variations in arrival times is PSR B1828-11 (Stairs, Lyne & 
Shemar 2000; Stairs et. al. 2003). Fourier analysis of these 
variations reveals harmonically related periodicities at ap- 
proximately 1000, 500 and 250 days (Stairs, Lyne & She- 
mar 2000; Fig. with the latter two somewhat more pro- 
nounced than the first. The length of the timescale of these 
variations implies that they are probably not of magneto- 
spheric origin, since the natural timescale in the magneto- 
sphere is of the order of the spin period, which in this case 
is P* = 0.405 sec. Even the ExB drift of subpulses (e.g. 
Ruderman & Sutherland 1975) does not exceed ~ 10P*. 
(However, Ruderman [2001] has suggested the possibility of 
drifts with periods of the order of a year.) As of the time of 
writing, there are no quantitative models for the data based 
solely on magnetospheric effects, but there are successful 
models based on precession (e.g. Jones & Andersson 2001, 
Link & Epstein 2001, and Wasserman 2003). Link & Ep- 
stein (2001) previously modelled the timing residuals from 
this pulsar in terms of precession of an axisymmetric, oblate 
rotating rigid body slowing down according to the vacuum 
magnetic dipole radiation formula. They found that the ob- 
servations could be accounted for in this model provided 
that the underlying pulsar is nearly an orthogonal rotator 
(magnetic obliquity \ ~ 89° to the body's symmetry axis) 
and nearly aligned angular momentum (wobble angle 9 m 3° 
between angular velocity and symmetry axis). These are in 
accordance with the conclusions reached by Jones & Ander- 
sson (2001). Although the precession amplitude is small, it 
may suffice to unpin superfluid vortex lines (Link & Cut- 
ler 2002), thus avoiding a potential impediment to preces- 
sion: pinning was shown to shorten the precession period to 
about 100 spin periods, and precession itself is dissipated 
over a timescale of 100-10,000 precession periods (Shaham 
1977, 1986; Sedrakian, Wasserman & Cordes 1999). Wasser- 
man (2003) argued that the data could also be accounted 
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Figure 1. Definition of various angles: the wobble angle (6) is 
the angle between the angular momentum (L) and the body z 
axis (which is chosen as the principal axis corresponding to the 
largest moment of inertia, J3); the beam swing angle (1?) is the 
angle between the angular momentum and the magnetic axis (b); 
X is the polar angle of the magnetic axis in the body frame. Note 
that the angles may not be coplanar. For an axisymmetric body 
the wobble angle remains constant; for a triaxial body it varies 
with time (see Appendix B). 



for if the underlying neutron star has either a type II su- 
perconductor or a strong toroidal magnetic field in its core. 
In these models, the angle between the angular velocity and 
symmetry axis of the star could be larger than found by 
Link & Epstein (2001), and the star does not have to be a 
nearly orthogonal rotator; spindown variations were found 
to dominate the timing residuals in this case as well. Link 
(2003) showed that the standard picture of the core of type 
II superconducting protons coexisting with superfluid neu- 
trons is inconsistent with long-period precession; pinning of 
the neutron vortices to the proton flux tubes makes the pre- 
cession frequency comparable to the rotation frequency of 
the star, a factor of 10 s too fast. Possible implications in- 
clude a normal core (both neutrons and protons), superfluid 
neutrons and normal protons, normal neutrons and super- 
fluid protons (type I or type II), or superfluid neutrons and 
type I protons. Sedrakian (2005) studied the last possibil- 
ity. He calculated the drag on neutron vortices moving in 
a type I superconducting core, and found that the drag is 
sufficiently small (that is, the vortices are sufficiently mobile 
with respect to the protons) that long-period precession is 
indeed possible in this scenario. Irrespective of the details, 
magnetic stresses in excess of the relatively weak ones that 
would arise from the pulsar's apparent dipole field strength, 
together with crustal stresses, would render the neutron star 
effectively triaxial in shape (Cutler 2002; Wasserman 2003; 
Cutler, Ushomirsky & Link 2003). 

Link & Epstein (2001) and Wasserman (2003) gave two 
alternative models that interpret the timing of PSR B1828- 
11 as precession. These models fit the data well, thus pro- 
viding strong evidence that the observed timing variations 
do indeed represent precession. These two models, however, 
are special cases. The purpose of this paper is to do a thor- 
ough search of the parameter space to see what we can learn 
about the properties of the spindown torque and the stellar 
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figure. To this end, we analyze the period residuals from this 
pulsar in terms of a simple model in which the rotating neu- 
tron star is assumed to be a triaxial rigid body. Obviously, 
precise axisymmetry is a special case, and we do not expect 
it to hold generally, particularly if the crust of the star is not 
in a relaxed state (e.g. Cutler, Ushomirsky & Link 2003), or 
has substantial internal magnetic stresses that may not be 
axisymmetric to begin with. Thus, one of our goals is to see 
what the data from PSR B1828-11 reveal about the shape 
of the neutron star crust. 

In this paper, we model a precessing neutron star as 
a single (rigid) body, rotating uniformly. Realistic neutron 
star modelling should take account of at least two different 
components - its solid crust, and (super)fluid core. Bondi & 
Gold (1955) considered the precession of a body consisting 
of a solid crust coupled frictionally to a fluid core. Their 
work showed that the long-term precession of the composite 
system depends on the timescale t cc on which the crust and 
core couple to one another. If Q.t cc <C 1 then the crust and 
core are very tightly coupled to one another on timescales 
smaller than a rotation period, and the moment of inertia 
tensor relevant to precession is that of the entire system, 
crust plus core. In this case, precession damps out slowly, on 
a timescale ~ (fltcc)^ 1 precession periods. If r2£ cc ^> 1 the 
crust and core only couple on timescales long compared to a 
rotation period, and the moment of inertia tensor relevant to 
precession is that of the crust alone. In this case, precession 
also damps slowly, with a characteristic decay timescale ~ 
Qt cc precession periods. Estimates of the crust-core coupling 
timescale vary but the consensus is that the coupling is weak, 
with Qt cc ~ 10 2 - 10 4 > 1 (e.g. Alpar, Langer & Sauls 1984, 
Alpar & Sauls 1988, Sedrakian & Sedrakian 1995), so the 
precession dynamics are governed by the moment of inertia 
tensor of the crust alone. However, it also turns out that 
as long as the crust and core couple on a timescale short 
compared with a precession period, but still long compared 
to the spin period, the relevant moment of inertia for all 
spindown effects, including those that vary periodically over 
a precession period, is the total stellar moment of inertia 
(Akgun, Link & Wasserman 2005, in preparation); this is the 
appropriate regime as long as Slt cc <C 10 8 , which appears to 
be the case. Thus, our one component model for precession 
is justified, apart from slow decay of the precession, which 
we neglect. 

Another issue is that the neutron star crust is not per- 
fectly rigid, but has a finite shear modulus. For a biaxial pre- 
cessing star, the crust must be strained in order for the star 
to precess with a period of order a year (Cutler, Ushomirsky 
& Link 2003). In addition, the strain field will vary with time 
as the star precesses, making the star's moment of inertia 
tensor time-dependent in the body frame. For simplicity, wc 
neglect these effects, and assume that the rotation of an im- 
perfectly rigid, triaxial star is well-described by the Euler 
equations for a rigid body, but with a moment of inertia 
tensor that is rescaled to account for the finite shear modu- 
lus. 

Since = 0.405 seconds and i sd w 10 5 years for PSR 
B 1828-11, F sd ~ 10 3 , and the spindown contribution to 
the precession-induced timing residuals is particularly im- 
portant. As a result, we may also hope to use these data to 
probe the value of a, that is, to probe the angular depen- 
dence of the spindown torque. While it is common to as- 



sume that a = 1 for rough analysis, on theoretical grounds 
we should not expect this to be true, for even an aligned 
rotator surrounded by a magnetosphere radiates energy, a 
process whose source is ultimately the rotational energy of 
the star, resulting in spindown at a rate presumably not 
much different from its luminosity divided by its rotational 
frequency. One of our chief findings is evidence that the ex- 
ternal torque that spins down a pulsar does indeed possess 
at least some angular dependence. 

Our analysis uses the same segment of the data 2 that 
was the basis of Link & Epstein (2001); this facilitates direct 
comparisons between the results of the two studies. We focus 
on the period residuals because we can derive an analytic for- 
mula for them in terms of elliptic functions. (An analogous 
formula for an oblate axisymmetric rotator has been derived 
previously by Bisnovatyi-Kogan, Mersov & Sheffer 1990, and 
Bisnovatyi-Kogan & Kahabka 1993, and has been applied to 
the 35-day cycle of Her X-l.) Because the underlying triax- 
ial model involves numerous parameters, using an analytic 
formula speeds up the computation considerably, which is a 
distinct advantage. By contrast, direct analysis of the timing 
residuals would require numerical integration of the model 
equations, a distinct disadvantage. Thus, for computational 
convenience, we analyze the period residuals rather than the 
arrival time residuals. However, a straightforward analysis 
of the period residuals using their tabulated uncertainties 
yielded very large values of \ 2 (~ a f ew x 10 3 ) under the 
assumption that the residuals are due solely to precession. 
This indicated to us that there is extra noise in the period 
residuals, either because their estimated uncertainties are 
too small (which we regard as unlikely to account for all the 
noise) or because there is a physical source of period noise 
that smears out the smooth contribution from precession 
systematically. In order to account for this extra noise sim- 
ply, we multiplied the tabulated uncertainties by a (single) 
factor F, and then marginalized over F to obtain posterior 
distributions of the (more interesting) parameters of the pre- 
cession model. This method - Student's t-test - represents 
a computational realization of "chi-by-eye" for data whose 
uncertainties may only be known incompletely. Details are 
given in Appendix D. 

Section 2 contains basic features of our model; further 
details may be found in Appendices A, B and C. Section 
3 contains results and implications of our analysis; statis- 
tical details, including the "chi by eye" method mentioned 
above, are found in Appendix D. Section 4 is a short di- 
gression on the pulse shape of PSR B1828-11, which is seen 
to vary systematically with precession phase (Stairs et al. 
2000). Although we do not use this information in our sta- 
tistical analysis, precession samples different regions in a 
pulsar's radio-emitting region, and offers the possibility of 
mapping out its shape (as has been done for PSR 1913+16, 
which exhibits geodetic precession, by Weisberg, Romani & 
Taylor [1989], and previously by Link & Epstein [2001] for 
PSR B1828-11). 



2 We thank I. H. Stairs, A. G. Lyne and S. L. Shcmar for gener- 
ously sharing their timing residual data with us. 



4 T. Akgiin, B. Link and I. Wasserman 




i i i i i ± i 

49,500 50,000 50,500 51,000 

Modified Julian Date 



Figure 2. Time of arrival residuals, period residuals, and shape 
parameter for PSR B1828-11 (courtesy of I. H. Stairs, A. G. Lyne 
and S. L. Shemar). The shape parameter is defined in terms of 
the weight of the narrow (A^) and wide (A^) standard pulse 
profiles that are present at every epoch, S = Apf/(Aj^ + Ay/), 
so that S ~ 1 for narrow pulses, and S ~ for wide ones (Stairs, 
Lyne & Shemar 2000; Stairs et. al. 2003). The solid line uses a 
cubic spline to connect the points for AP. 



Table 1. Definitions of important parameters. 



parameter physical meaning 

, polar and azimuthal angles of the 

magnetic axis in the body frame 

e 2 measures the degree of triaxiality 

determines the components of the 
A angular momentum and is 

related to the wobble angle 

determines the strength of the angular 
part of the spindown torque 



2 MODELS 

Here we review the models that we use briefly, highlighting 
some of the more important parameters. The derivations for 
the period residuals are lengthy and are left for the Appen- 
dices; we present the geometric model in Appendix A, and 
the spindown model is derived in Appendix C. In what is 
next, we will follow the notation used in these Appendices. 
The parameters of our models are listed in Tabled 

2.1 Geometric Model 

What we refer to as the geometric model is the effect of 
triaxiality alone (i.e. torque-free precession). In this case, 
Euler's equation for the angular momentum can be solved 
analytically in terms of Jacobian elliptic functions (Landau 
& Lifshitz, 1976). As we show in Appendix A, the period 
residuals are then found to be of the form AP ge /P* ~ w v f n 
where P* is the rotation period at a fiducial epoch, vj p is 
a dimensionless quantity of order P*/P P ~ e, and f n is a 
complicated combination of the elliptic functions. Because 
of the inherent form of f n the amplitude of the residuals is 
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Figure 3. Sample models for the period residuals. The data 
points are shown for comparison. The solid line is for the full 
model (both effects included), the dashed line is for the geomet- 
ric model alone, and the dotted line is for the axisymmetric model. 
The parameters for all models shown here are quoted in Table U\ 



Table 2. List of parameters for the models shown in Fig. l3l 



parameter 


geometric 


axisymmetric 


full 


X 


74.0° 


88.6° 


71.8° 


4> 


12.0° 


0° 


0° 


e 2 


3912 





2135 


A 


0.012 


0.0437 


0.00325 


a 





1 


0.983 




18.5° 


5° 


0.44° 



not trivial to predict in general, and they can exhibit very 
rich behavior. 

We denote the principal moments of inertia by Ii\ 
the associated axes serve as the basis for the rotating 
(body) frame. We then define the following parameters: 
e = (Jz — which measures the deviation from spheric- 

ity; e 2 = [h{h-h)]/[h{h-h)} which measures the degree 
of triaxiality; k 2 which is the parameter of the Jacobian ellip- 
tic functions, and depends on the angular momentum and 
the moments of inertia; and A which determines the com- 
ponents of the angular momentum (and would be simply 
A = L1/L3 for an axisymmetric star, but is slightly different 
in the more general case, see Eq. IA4I 1. The last three are 
not independent: k = eA. Note that k 2 does not depend on 
e, but only on e 2 and A. 

Note that for e 2 = the body is oblate and axisymmet- 
ric, and A = when there is no precession. In both cases 
k — 0, and the Jacobian elliptic functions reduce to the reg- 
ular trigonometric functions. At k = 1 they become hyper- 
bolic functions, and the angular momentum exponentially 
aligns with the principal axis corresponding to the interme- 
diate moment of inertia, Ii- Between these two extremes, 
the precession of the angular momentum takes place along 
the intersection of the sphere defined by the conservation of 
angular momentum (L 2 = LiLi), and the ellipsoid defined 
by the conservation of energy (E — L 2 /2/;). The result- 
ing shape of the trajectories is the Binet ellipsoid (Landau 
& Lifshitz, 1976). In the limit e 2 — > 00 the body becomes 
prolate axisymmetric. 
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2.2 Spindown Model 

The rotation of an isolated neutron star is not torque-free, 
but slows down, resulting in a gradual increase in the rota- 
tion period. It is thought that because of the rotating mag- 
netic field, angular momentum is lost to radiation. In the 
simple model of a rotating dipole in a vacuum, the pulses 
are emitted at the poles of the magnetic field and the torque 
has the form N ~ —N [Cl — (f2 ■ b)b], where ft is the in- 
stantaneous rotation axis, and b is the pulse (and magnetic) 
axis. 

This is a very crude model. The magnetic field may have 
non-dipolar components of considerable amplitude, which 
we will not consider here. The pulsar is also not in a per- 
fect vacuum, but is surrounded by a plasma-filled magne- 
tosphere (Goldreich & Julian, 1969). The vacuum torque 
vanishes when the angular velocity and the magnetic axis 
are aligned, while the presence of a magnetosphere would 
require a loss of angular momentum, no matter what the 
orientation is. Therefore, the vacuum dipole torque should 
give only an incomplete description at best. We adopt a gen- 
eral spindown torque of the form N s d = — N [fl — a(tl -b)b], 
where we introduce an additional parameter a. Loss of an- 
gular momentum mandates that a ^ 1. It should also be 
positive, or we would have an angle dependent torque that 
is opposite in sign to the dipole contribution. The vacuum 
case is retrieved by setting a = 1. The case of a = corre- 
sponds to an external torque with no angular dependence, 
which would not produce periodic time of arrival residuals. 

The torque-modified Euler's equation can be solved ap- 
proximately for small e, and a second contribution to the 
period residuals arises due to the torque, AP s d/P± ~ T a dA£. 
We will refer to this as the spindown model; and the sum of 
both geometric and spindown contributions will be referred 
to as the full model. Here, T s d = IzNo/vjpL 2 ~ P P /t B d ~ 
10 -5 and is determined by the spindown properties of the 
neutron star (in particular, the period derivative, P or the 
characteristic age, t s d); Al is another complicated function 
of Jacobian elliptic functions and Legendre integrals. For an 
axisymmetric star, Al = ai sin ui p t + a% sin 2uo p t, where ui p is 
now the precession frequency, and ai are some coefficients; 
the axisymmetric model thus has two harmonically related 
components. We take the two peaks in the spectrum of the 
period residuals of PSR B1828-11 with periods of ~ 500 and 
~ 250 days to be the most significant. From the form of the 
axisymmetric model, we thus conclude that the precession 
period must be ~ 500 days. It can also be shown that for an 
axisymmetric star, in the region of interest, the geometric 
contribution is quite negligible compared to the spindown 
contribution for a ~ 1 (see Appendix C). 

If the 1000-day period represented the precession pe- 
riod, we would expect to see variations of the pulsar beam 
width at the same period. While such changes were reported 
by Stairs, Lyne & Shemar (2000), subsequent more detailed 
analysis has not shown a 1000-day period in the beam width 
data (Parry et al. 2005). We thus assume that the precession 
period is ~ 500 days, and attribute the 1000-day component 
in the timing data to something unrelated to precession, such 
as timing noise. We note that our model cannot provide sat- 
isfactory fits to the data if the precession period is ~ 1000 
days. 



2.3 Constraints and Statistical Analysis 

The orientation of the angular momentum, L is fixed in the 
inertial frame. This is still true even in the presence of the 
spindown torque, if e is sufficiently small (see Appendix C). 
Then the requirement that the pulse beam, which we assume 
to be centered along the magnetic axis, b never precesses 
entirely out of our line of sight means that the angle between 
L and b should not vary by more than the angular width of 
the pulse itself. We will refer to the angle between L and b 
as the beam swing angle and denote it by If the angular 
radius of the pulse is p, then the above constraint can be 
expressed as Ai? = i? maa; — t? m in ^ 2p; in general, we will 
require the beam swing variation to be less than some value 
Afimax- We also will define the wobble angle, 9 as the angle 
between the angular momentum and the body z axis (see 
Fig- and Appendix B). 

The duty cycle allows us to estimate the angular extent 
of the pulse, and for PSR B1828-11 this varies between 5° 
to 7°. For a circular pulse, this implies that the beam swing 
angle cannot be varying by more than a few degrees. Larger 
variations would require a more elongated pulse shape. Yet, 
at this time, there is not enough evidence to elaborate more 
on this. In particular, polarization data might be quite useful 
to determine the extent of the pulse. Stairs et al. (2000) also 
report periodic variations in the average pulse shape. We 
offer a possible explanation in a following section. 

Another restriction may be that PSR B1828-11 does 
not have an interpulse. That can further restrict the rela- 
tive orientation of the angular momentum and the magnetic 
axis. However, due to uncertainties in the structure of the 
magnetic field, it is not clear that an interpulse will neces- 
sarily appear. Therefore, we do not impose this restriction. 
The observer's location is an additional parameter, and can 
be independently fixed. 

We apply the two models - geometric (a = 0) and full 
(a / 0) - under the given constraints to PSR B1828-11, us- 
ing a Bayesian approach to obtain probability distribution 
functions (pdfs) for individual parameters. We assume spe- 
cific priors in the full multi-dimensional parameter space (to 
be discussed next), but the effective priors exhibited in the 
projected (marginalized) 1-D posterior pdfs shown in the fig- 
ures below are integrals over the constraints. Because there 
appears to be systematic noise in the period residuals larger 
than their tabulated uncertainties, we scale the latter and 
then marginalize over the scaling factor as detailed in Ap- 
pendix D. Once the likelihood is determined, the individual 
pdfs are obtained through integration over the remaining 
parameters and normalization. 

Let {pk} denote the set of the n parameters. Then the 
likelihood, £({p k }) and the volume of integration, V{{pk}) 
are functions of this set. The latter also depends on the beam 
swing angle constraint, which itself is a function of a subset 
of the parameters. The priors, gi{pi) are functions only of 
the single parameter they refer to. Then the projected 1- 
D posterior pdf for the i-th parameter can be expressed as 
an integral of the likelihood over the remaining parameters, 
over the volume defined by the constraints, 

p n 

fi({Ph})= gi(pi)£({Pk})T\gi(pj)dpj . (l) 
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Figure 4. Pdfs for beam swing angle variation less than 1°. In 
Figs. IHSl the prior for \ i s A a t over d(cosx); all other priors are 
flat. The cutoff for e 2 is at 2, and A is confined to be below 0.2. x 
is in degrees; a, e 2 and A are dimensionless. The percentages listed 
above the plots are the probabilities for the full model enclosed 
in the corresponding ranges. 

Similarly, the projected 1-D prior can be expressed as, 

p n 

hi{{Pk})= / g l {j>i)WgjijP]) d P] ■ (2) 

Jv{{ Pk }) m 

It is these two quantities (/; and hi) that are plotted in Figs. 
1 111)1 Note that, if the volume of integration had not depended 
on the constraints, then we would simply have hi = g%. 

Within the context of our precession model we can use 
the pdf to compare how well different sets of model pa- 
rameters fit the data. However, we cannot assess the extent 
to which the data demand explanation in terms of preces- 
sion, as opposed to some other, completely different physical 
model. Any model for the data will lead to a pdf with local 
maxima at certain values of the parameters of the model, 
and we can assess the relative significance of these peaks to 
quantify the extent to which the model parameters are de- 
termined by the data. Whether or not another model that is 
just as well-motivated physically as our precession model can 
fit the data better is outside the scope of our analysis. Given 
a competitor model - of which we are unaware - Bayesian 
methods could be used for making model comparisons. 



3 RESULTS AND DISCUSSION 

The physical parameters that determine the form of the 
residuals are the two angles that specify the orientation of 
the magnetic axis in the body frame (the polar angle, \ an d 
the azimuthal angle, tp; see Fig. IA1L any two of e 2 , A and 
k 2 ; and a. There is also a r (measured in precession cycles) 
that determines the initial phase. Thus, the total number 
of parameters is six. \ varies between and tt/2, and its 
prior is taken to be flat over cosx; an d <j> varies between 
and 2n, and has a flat prior. In other words, we assume that 
orientations of the magnetic axis are equally likely over all 
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Figure 5. Pdfs for beam swing angle variation less than 3°, for 
the same set of priors as in Fig. El 
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Figure 6. Pdfs for beam swing angle variation less than 5°, for 
the same set of priors as in Fig. |3] 

solid angles. Priors for t and a are flat between and 1. 
On the other hand, e 2 and A can have any positive values, 
as long as the beam swing angle is constrained and k 2 < 1; 
therefore, we have to introduce cutoffs in their priors. A is 
related to the wobble angle, and due to the beam swing an- 
gle constraint it cannot be too large; we take A ^ 0.2 with 
a flat prior. 

The situation is slightly more complicated for e 2 . The 
crust of a neutron star (which in our model is the only com- 
ponent since we do not consider the liquid interior) can relax 
only through shearing motions as it spins down and so must 
be triaxial (Link, Franco & Epstein 1998; Franco, Link & 
Epstein 2000). Adding the magnetic stresses, which result 
from the multi-polar field near the surface would produce a 
very complicated figure. It is, therefore, quite unlikely that 
the star is oblate axisymmetric (e 2 = 0) or prolate axisym- 
metric (e 2 — > oo) to a very high precision. On theoretical 
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Figure 7. Pdfs for beam swing angle variation less than 1°. In 
Figs. 17151 the prior for \ is fiat over cosx; the prior for e 2 is flat 
over ln(l + e 2 ); all other priors are flat. The cutoff for e 2 is at 
4000, and A is confined to be below 0.2. \ i s m degrees; a, e 2 and 
A are dimensionless. e 2 is plotted on a log-log scale to reveal more 
detail. The percentages listed above the plots are the probabilities 
for the full model enclosed in the corresponding ranges. The prior 
probabilities are given for e 2 in parentheses for comparison with 
the posterior probabilities. 



grounds one might expect e 2 to be close to unity. Thus, we 
first take e 2 ^ 2 with a flat prior (Figs. I HOP . However, we 
find that within this range the pdf for e 2 is not confined. Be- 
cause of that, we also consider a second case where we allow 
for larger values of e 2 (Figs. 17191 . The spindown model we 
use is derived under the assumption that e 2 is not exceed- 
ingly large (see Appendix C). Therefore, we take e 2 ^ 4000, 
which seems to encompass the regions of interest, without 
violating our assumptions. Since the volume of integration 
is considerably larger at large values of e 2 than at small 
values, taking a flat pdf over e 2 in this case would greatly 
suppress the importance of small e 2 . Therefore, we need to 
incorporate a prior that is fair for both regimes: we take a 
flat prior over ln(l + e 2 ), i.e. the prior for e 2 is l/(l + e 2 ). For 
both e 2 < 2 and e 2 < 4000, we calculate pdfs for three dif- 
ferent values of the maximum beam swing angle constraint: 
A3 max = 1°, 3° and 5°. 

In Fig. [3] we show the data that we use, together with 
some sample models. The parameters for these models are 
listed in Table 1. The best fit that we find is a purely ge- 
ometrical model, which has a very large beam swing angle 
that is, in fact, outside our prior range (which was relaxed 
for determining an unconstrained, global "favorite" model). 
The axisymmetric model given here is similar to that of Link 
& Epstein (2001), except that the beam swing angle is con- 
strained to be below 5°; in fact, as we show in Appendix 
C, any axisymmetric model is assured to yield quite sim- 
ilar results even when we introduce the additional torque 
parameter a. 

Because of the large number of parameters, numerical 
integration for the pdfs is quite time consuming. The fig- 
ures presented here typically have a resolution of about 100 
points per parameter or less. This means that fine structure 



Figure 8. Pdfs for beam swing angle variation less than 3°, for 
the same set of priors as in Fig. 171 
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Figure 9. Pdfs for beam swing angle variation less than 5°, for 
the same set of priors as in Fig. |7J 



in the pdfs may have been missed. Nevertheless, the most 
notable structures in the pdfs are expected to remain. 

In Figs. l4l9l we show the projected 1-D pdfs for each pa- 
rameter, computed by integrating the multidimensional pdf 
over all other parameters. The prior is also a function of the 
entire set of parameters and is not separable for all except a. 
We define the 1-D prior for a given parameter by integrating 
over the rest. A comparison with the full 1-D pdf illustrates 
the importance of the period residuals in determining the 
pdf. Keep in mind that both the prior and the posterior 
pdfs also include the beam swing angle constraint. In these 
figures, the dotted lines are for the prior pdfs, the dashed 
lines are for the geometric model alone (Eq. HA27I ). and the 
solid lines are for the geometric model and the spindown 
model (Eq. 1C36H both combined. 

We now discuss some of the main characteristics and 
implications of our analysis. 
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Figure 10. Multivariate pdf surfaces for the full model at constant 
are integrated out, and normalization is carried out over all surfaces 

The torque parameter a: For e 2 ^ 2 Figs. l4lb1 show con- 
siderable probability over the whole range of acceptable val- 
ues, with a peak at low a that becomes more prominent 
as Ai? increases. Another lower and wider peak appears for 
Ai? maa , = 5° at larger values of a (Fig. |5J. Nevertheless, 
neither of the peaks is highly significant, because they do 
not contain most of the probability. Therefore, we conclude 
that the data do not constrain a strongly, and it can be 
quite different from the vacuum spindown value a = 1. As 
larger values of e 2 are permitted, Figs. 1 7191 show a peak at 
a-tl, but with a large tail extending over most of the pa- 
rameter space. With increasing values of A#, lower a values 
become likelier, but most of the probability (> 90%) still 
lies at a 0.25. Thus, the value of a is not well-determined, 
but there is evidence for an angle-dependent torque. The 
parameter a is truly a measure of the angle dependence of 
the spin-down torque; it does not depend on the geometric 
effect at all. 

The magnetic inclination \ : The axisymmetric model 
discussed by Link & Epstein (2001) requires \ to be ex- 
tremely close to 90°. As discussed in Appendix C, this is 
true even when we allow a 7^ 1. For a triaxial model, we 
find the range of acceptable \ values to be much larger. 
The geometric model has a peak at small \y which moves 
on to higher values of \ an(1 broadens with increasing beam 




e 2 as functions of a and \ for Ai? < 3° . The remaining parameters 
. The ripples are artifacts of integration. 

swing angle. This trend is still seen at e 2 ^ 2 when a / 
is turned on. The spindown produces a narrow but strong 
peak in the vicinity of 90° , which becomes more pronounced 
as A$ is allowed to be bigger. This peak corresponds to the 
axisymmetric case, and implies that it requires larger At? 
values; in fact, the model discussed by Link & Epstein has 
Atf ~ 6.4°. For e 2 < 4000 (Figs. EE), the peak at large x re- 
mains apparent, though now it is quite broad. The inclusion 
of points beyond e 2 ~ 2 seems to favor a more important 
spindown contribution, and the geometric effect is further 
suppressed. There is also a small cusp that appears in the 
pdf for Ai? < 5°, at x very near 90°, corresponding to the 
axisymmetric case. Yet, this peak is quite narrow, and the 
vast majority of the probability lies outside of it. 

The triaxiality parameter e 2 : For e 2 ^ 2, the pdf looks 
quite similar to the prior, implying that the data do not dif- 
ferentiate among values of e 2 (Figs. 14161 . There is a narrow 
sharp peak at e 2 = for the full model, corresponding to 
the oblate axisymmetric case, but the probability enclosed 
within the peak is very small. Note that, at the same time, 
the pdf for the geometric effect alone almost vanishes, i.e. 
for the axisymmetric case, spindown is essential. When we 
allow e 2 > 2, another peak appears in the pdf (Figs. 17191 . 
At these values of e 2 the star is once again almost axisym- 
metric, except that now 7i < I2 — I3, i.e. the star is prolate. 
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Qualitatively e 2 = 1 separates oblate and prolate shapes. 
Then, comparing the probabilities enclosed in the two re- 
gions, e 2 < 1 and e 2 > 1 , we find that the prolate case con- 
tains, by far, most of the probability, even though we have 
adopted a prior which somewhat disfavors large e 2 values. 
Note that as the beam swing angle constraint is relaxed, the 
probability becomes quite evenly distributed over a wide re- 
gion in e 2 . This leads to the conclusion that there are many 
triaxial models with a wide range of e 2 that can fit the data. 
In other words, the data do not discriminate among values 
of e 2 , especially when Ai9 is relatively large. 

The wobble parameter A: For both e 2 ^ 2 and e 2 
4000, the pdf is contained in a region which seems to be 
mostly confined by the beam swing angle. The constraint 
results in a dramatic cutoff at the high end of A. Beyond 
that point, we cannot find a value of e 2 for which A# will 
be smaller than the constraint. Within that region, there is 
considerable probability distributed over the whole range of 
A; both models follow the same trend. We can conclude that 
the oblate case favors somewhat larger values of A, which 
produces a peak that is partially visible for A-d ma x = 5° 
(Fig.®. The prolate case, on the other hand, favors smaller 
values of A, resulting in a second peak, which appears when 
we allow e 2 to be large (Fig. EJ, but is absent when e 2 is 
confined to low values (Fig. 0. 

Fig. HOI shows multivariate pdfs for the full model, plot- 
ted as surfaces at constant e 2 , as functions of a and X- The 
remaining parameters (0, A and t ) are integrated out, and 
the pdfs are normalized over all surfaces. (The special case 
e 2 = can be done with considerably higher resolution, be- 
cause of the freedom of choice of <f>.) The ripples that are 
present are artifacts of the integration and subside as the 
resolution is improved. The amplitude of the ripples serves 
as an implicit way of evaluating the significance of the peaks 
in the pdf. Relatively large amplitude implies that there are 
no significant peaks in the pdf, meaning that no points or 
regions in the allowed parameter space are favored strongly. 

As Ai9 mal is increased, two regions acquire prominence: 
a very narrow ridge at large Xi which extends over a wide 
region in a and includes the axisymmetric case, and a smaller 
peak at small \ an d small a. Keep in mind that this second 
peak is further discriminated against by the prior, which is 
flat over cosx, i.e. there is a factor of sinx that also enters 
the pdf. Consequently, when the one-dimensional pdfs are 
calculated, the second peak is considerably suppressed. 



4 SHAPE OF THE PULSE 

The pulse profile of PSR B1828-11 alternates between two 
different modes, one narrow and the other broad, and has 
Fourier power at both 250 and 500 days (Stairs, Lyne & She- 
mar 2000, Stairs et al. 2003). Stairs et al. (2003) describe 
how the pulse profile is determined. During a particular ob- 
serving session, 16 pulse averages may be either broad or 
narrow, with a shape parameter S defined to be the fraction 
of the mean pulse shape for that session attributed to the 
narrow component. The shape parameter 5* varies system- 
atically between « (all wide) and « 1 (all narrow) over 
the precession cycle, with a strong Fourier component at the 
"first harmonic" 1/250 days of the "fundamental" precession 
frequency 1/500 days (Stairs, Lyne & Shemar 2000, Stairs et 




Figure 11. Schematic of a pulse consisting of a bright core, and 
surrounded by smaller fainter conal blobs. Core emission is as- 
sured to be stationary, whereas conal emission could vary as the 
emitting blobs circulate about the beam axis. 



al. 2003). Link & Epstein (2001) suggested that the emission 
beam of the pulsar must have an hourglass shape in order for 
S to exhibit substantial variability on the 250 day timescale. 
(A similar elongated shape was inferred from studies of the 
geodetic precession of PSR B1913+16 by Weisberg & Tay- 
lor 2002.) However, they did not address the issue of mode 
switching during individual observing sessions at all. Here 
we present an alternative viewpoint centered around model- 
ing profile mode switching within individual observing ses- 
sions, and argue that it may be able to produce some aspects 
of the required harmonic structure shown in Fig. 

The basic geometrical picture is shown in Fig. 1111 We 
attribute the narrow component of the pulse to core emis- 
sion centered around the beam axis. The broader profile is a 
superposition of core and conal emission. Thus, in the par- 
lance of Rankin (1990, 1993), the pulsar alternates between 
presenting a core single (St) and triple (T) pulse profile. The 
relatively young spindown age of PSR B1828-11 (about 0.11 
Myr) and its large value of B12/P 2 are consistent with this 
categorization. However, the apparent pulse width in the 
narrow state appears to be anomalous: Rankin (1990) finds 
a FWHM pulse width W COIC = 2.45° P~ 1/2 / sin X , where \ 
is the angle between a pulsar's spin and magnetic axes. For 
PSR B1828-11, the FWHM of the narrow state is about 
2.3°, as opposed to Wcore ~ 3.85°/ sinx from Rankin's 
(1990) formula. We note that the bounding relationship 
W CO ic ^ 2.45°/P _1//2 was derived from a set of "interpul- 
sars" thought to be nearly orthogonal rotators, so we would 
have expected Rankin's (1990) formula to work especially 
well if x is near 90°. as was suggested by Link & Epstein 
(2001); on the other hand, the discrepancy is also small- 
est for x ~ 90° , which may be circumstantial evidence that 
PSR B 1828-11 really is nearly an orthogonal rotator. (More- 
over, the frequency dependence of the core width is relatively 
weak at high frequencies, so that the fact that the Rankin 
[1990] formula is for 1 GHz emission, whereas the Stairs et al. 
[2003] observations were at 1.6 or 1.7 GHz is not responsible 
for the discrepancy.) We note that there are other excep- 
tions to Rankin's (1990) bound, but not many (see e.g. Fig. 
23b in Graham-Smith [2003], adapted from Gould [1994]); 
given uncertainties in X there may be other pulsars with 
W core > 2.45 / VP but also Wcorc < 2.45°/ VP sin X - 

If pulsar core emission intensity were Gaussian, we 
would expect an observed intensity of the form, 
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-Tcore(O) exp 



l_ 

2p\ 



2pl 



(3) 



where is the impact parameter of the observer's line of 
sight relative to the beam axis, <j> is pulse phase (centered 
on epoch of closest passage relative to the axis), and a = 
X + P is the angle between the line of sight to the observer 
and the stellar spin axis. Here pi and p2 define the extent 
and the shape of the beam, which would be elliptical when 
they are not equal. This formula assumes that /3 <C \, and 
that emission is strongly beamed along magnetic field lines, 
but does not presume that the emission pattern is circularly 
symmetric with respect to the beam axis. The two directions 
1 and 2 are relative to a coordinate system in which £3 = b 
coincides with the magnetic moment of the star, which is 
assumed to be the beam axis, £2 = Lxb/|Lxb|, and ei = 
(b-Lb — L)/|Lxb|. For a Gaussian beam, Eq. shows that 
the core component width is independent of the impact angle 
(5, although the peak intensity is not (e.g. Rankin 1990). 
However, this is a unique property of a Gaussian profile. We 
can well imagine that the emission cuts off sharply (even 
discontinuously) for sufficiently large /3, in which case the 
core width could be narrower than normal. Because the peak 
core intensity would also be lower in such cases, it would 
be harder to detect, which may account for the rarity of 
exceptions to Rankin's (1990) bound. 

A sharper cutoff to the core emission beam would not 
only allow narrower core pulse profiles, but would also in- 
troduce p dependence into the width. As a simple example, 
suppose that the beam profile is, 
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i.e. still a self-similar function but with a sharper cutoff than 
a Gaussian profile. The peak intensity is at <j> = 0, where 
u = iimin = p 2 /pi; the FWHM is at phases ±<^i/2, where 
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(5) 



where the approximation is valid for small values of k, ir- 
respective of KUmi n . The cutoff becomes important once 
ftitmin ~ 1 i.e. for f3 > pi/y/it. The core width decreases 
with increasing f3, as does the peak intensity observed from 
the core. 

As we mentioned above, we ascribe the broader pulse 
profile state to a superposition of core and conal compo- 
nents. In keeping with the schematic Fig. 1111 we assume that 
the conal emitting pattern is patchy and, as we discuss fur- 
ther below, probably only stationary in the mean. Consider 
an individual conal emitting region (hereafter "blob") i; we 
assume that it is centered at (xj,x'f) = pi(cos<7i, sinai). 
The emission pattern of blob i may be anisotropic in a 
complicated fashion, with possible preferred directions not 
only along the ei,a axes, but also along and perpendicu- 
lar to Xj = (cos o"i, sinui). The observer sees an inten- 
sity that is a function of the two variables (3 — pi cos Oi 
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Figure 12. Acr/27r vs. /3 for various q and 8/p = 0.1. 
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Figure 13. "Shape parameter" vs. precession phase for a model 
with 8/p = 0.1, q = —0.8, N = 6, and impact parameter /3/p = 
0.95 + 0.15 sin <L. 



and — (<f) + 4>d) sin a — pismai, where the phase delay is 
<f>d = hi/cPi, w l°(/ii/340km), where hi is the height of 
blob i above the core emitting region. The peak value of 
the intensity of radiation seen from any given blob is only a 
function of \f3 — pi coscr;|, though, and we assume that blob 
i is detectable provided that 



|/3 — pi COS(Ti| < Si 



(6) 



where Si may depend on c,. Thus, the detectability of an 
individual blob varies through the precession cycle, and the 
probability of seeing any blobs at all also varies, thus affect- 
ing the observed beam width. 

The problem of modeling the detectability of a given 
blob can be quite complex. To illustrate, suppose that each 
blob is at p; = p, and has the same anisotropic shape. Sim- 
plify even more by assuming that the emission profiles of the 
blobs have a characteristic length 8 r along the (radial) direc- 
tion from the beam axis to its center, and a different length 
S t in the direction tangential to it. Then we can detect the 
blob if, 



(/3 — pcosCTi) ^ S t sin cr; + S r 



(7) 
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For convenience, we define 8 2 = |(<5 t 2 + <5 2 ) and q5 2 — \{8 2 — 
5 2 ,). Note that q may be positive or negative, and \q\ ^ 1; 
q — for a circular beam. Since we expect that in general 
5 <C p, and except for special values, S <C /3, we only expect 
blobs within a small range Aa(/3) to be visible. 

Fig. H2l shows the result of solving Eq. JJJ for the range 
of observable Aa/2n as a function of impact parameter (5. 
(The solutions were not extended beyond f3/p = 1+5^/1 — q, 
where Aa = 0.) The results exhibit complicated behav- 
ior even in this simple model. Given N conal blobs, the 
probability of seeing the broader pulse profile is large when 
2 x NAa/2-K > 1, and is small when NAa/n < 1. (The 
factor of two is because the observer's line of sight crosses 
the cone twice.) Thus, we may expect the shape parameter 
S to be small for impact parameters where Aa/2n is large, 
and vice-versa; during a precession cycle, both regimes may 
be sampled. 

Fig. ll3l shows an example of how the probability of see- 
ing only the narrow pulse would vary with precession phase 
in this model; this example captures the main features of 
the observed beam width variations shown in Fig. [5] For 
constructing the figure we adopted q = —0.8, S = O.lp, and 
assumed a sinusoidal variation of the impact parameter with 
precession phase: 

(3{<j> p ) =/3 + /3i sin <j> p (8) 

with Po = 0.95p and fi\ = 0.15p assumed for graphi- 
cal purposes. The "shape parameter" is taken to be S = 
(1 — Aa/iv) N i.e. the probability that no blobs are detected; 
Fig. 1131 assumes N = 6. Clearly, S varies periodically but 
not sinusoidally in this model, and also varies substantially 
in half a precession cycle. This distinctive "doubly periodic" 
variation is only seen if the observer's line of sight crosses 
near /3 — p. This is consistent with our earlier discussion 
of core widths if the core emission is still visible but start- 
ing to cut off at such impact angles. Presumably, the peak 
intensity of the core emission must also far exceed that of 
any conal blob for this model to be viable; there are some 
indications that conal emission becomes more prominent as 
pulsars age (e.g. Rankin 1990). Although the range of vari- 
ation of S in this example is smaller than in PSR B1828- 
11, extensions of the model, such as different assumptions 
about the conal emission (e.g. an hourglass-shaped cone as 
in Link & Epstein [2001], or a more complicated version of 
blob anisotropy) may possibly yield a better account of the 
data. 

If the conal blobs were stationary in the rotating frame, 
then the observer would see pulse profile variations as a func- 
tion of precession phase, but would not see any variations 
at a given precession phase. However, it is likely that the 
conal blobs are not at fixed positions but rather circulate 
around the cone in a rotating carousel (Deshpande & Rankin 
1999, 2001; see Rankin & Wright 2003 for a review). In this 
picture, which has empirical support (Deshpande & Rankin 
1999, 2001), conal emission is from beams that circulate with 
a frequency fid = /dfi* relative to a reference frame rotaing 
with the star; the circulation is probably the consequence 
of ExB drift (e.g. Ruderman & Sutherland 1975, Gil, Me- 
likidze & Geppert 2003, Wright 2003), and f d < 0.1 is a 
reasonable value. By contrast, the core emission is station- 
ary, and from a much lower altitude than the conal emission 
(possibly from near the polar cap). In this picture, during a 



given observing session the core component is always visible, 
but the conal component fluctuates as emitting blobs pop 
in and out of the observer's line of sight periodically. The 
probability that the observer sees conal emission at all varies 
systematically during the precession cycle, and is fixed dur- 
ing any observing session lasting a day or so (i.e. far less 
than the precession period) . 

If this model is correct at least in a broad outline, then 
the total beam swing during a precession cycle is sC 1 — 
3°, given expected core radii (Rankin 1993). Larger beam 
swings could be accommodated by a more complex model for 
the pulse shape (eg. the hourglass shape of Link & Epstein 
2001). 



5 CONCLUSION 

We find a wide range of triaxial models that may explain 
the period residuals of PSR B1828-11 in terms of preces- 
sion, even under the stringent constraints we have imposed. 
We find many fits that are as good as, or better than the ax- 
isymmetric model considered before (Link & Epstein 2001). 
In general, fits improve with larger beam swing angle varia- 
tions (A#), but if we assume that the pulse is confined to a 
region a few degrees in size, we have to rule them out. Pro- 
late and oblate axisymmetric models seem to be favored, 
especially for small A$, but that is not sufficient to rule out 
other triaxial models. Both the geometric and spindown ef- 
fects contribute to the fits. Oddly, if we relax our beam swing 
constraint completely, the data prefer a best fit that has 
a — (no spindown contribution), but Ai9 for that model is 
unreasonably large (Fig. |3] and Table 1), so it is merely an 
unphysical curiosity. 

In the oblate axisymmetric model, spindown is the dom- 
inant effect, but it requires parameters (in particular, a large 
X value) that could be expected to produce an interpulse, 
which is not seen in PSR B1828-11. If we were to impose the 
absence of an interpulse as a constraint, some of the models 
we have permitted in our analysis would be excluded, partic- 
ularly those at large x- Conceivably, the magnetic field and 
core beam structure of PSR B1828-11 are sufficiently com- 
plex that an interpulse would be absent even at x ~ * 90°. 
We note that our model for shape variations suggests that 
we are only viewing the outskirts of the core emission in the 
component we detect, which may enhance the probability 
that emission from the opposite pole is undetectable. Thus, 
we do not impose the absence of an interpulse as a constraint 
on our analysis. 

Our models do not require a = 1, so substantial de- 
viations from the vacuum spindown formula are allowed. In 
fact, rather small values of a are permitted for an oblate star 
(e 2 < 1). However, for a prolate star (e 2 > 1) we find that 
larger values are favored (a > 0.25), thus providing evidence 
for an angle- dependent torque. To our knowledge, our anal- 
ysis provides the first evidence that pulsars are spun down 
by a torque that depends on the angle between the magnetic 
moment and the instantaneous angular velocity. 

The magnetic obliquity, x i s no longer required to be 
extremely close to 90°, and we find A (which is related to 
the wobble angle) to be restricted mainly by the beam swing 
angle. Two peaks in e 2 are prominent, corresponding to the 
oblate axisymmetric (e 2 = 0) and prolate nearly axisymmet- 
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ric (large e 2 ) cases. These are especially evident for small 
beam swing angle variations. For larger beam swing varia- 
tions, the data do not discriminate among values of e 2 very 
much. 

In summary, our precession model fits the data equally 
well for a broad range of parameters. We cannot constrain 
the shape of the star, but we do find evidence for angle- 
dependent spin-down torque. Overall, the ability of our 
physically-motivated model to account for the principal fea- 
tures of PSR B 1828-1 l's timing data without special choices 
of the parameters reinforces the idea that the pulsar is pro- 
cessing. In particular, we have shown that the data can be 
fit without resorting to a nearly orthogonal rotator with 
a vacuum-like dipole torque as Link & Epstein (2001) did 
in their preliminary work. Though we cannot strongly con- 
strain the angular dependence of the spin-down torque with 
the present data, the potential remains for learning more 
about this important aspect of the neutron star magne- 
tosphere from future observations. The parameters of our 
model are not very tightly constrained. Two possible rea- 
sons are that we have only three cycles of data, and that 
the data have a large degree of intrinsic scatter that our 
simple model cannot account for, creating wide pdfs. It will 
be interesting to see if the parameters can be more tightly 
constrained as more data become available over the next 
decade. 

Our analysis does not employ the data on the shape 
parameter variations, because constructing a comprehensive 
mathematical description would require a reliable model for 
the pulsar beam. Although we do not possess such an accu- 
rate model, we offer an explanation using a compound pulse 
structure, with core and cone components. 

Precession has interesting implications for pulsar obser- 
vations, which so far have not been widely discussed. One 
immediate, and very obvious effect would be disappearing 
pulsars, i.e. pulsars that due to precession, would at some 
point leave the line of sight of the observer, but excluding 
other effects, would eventually come back. The timescale of 
such changes could be months to years. 



Deshpande A.A., Rankin J.M., 1999, ApJ, 524, 1008 
Deshpande A.A., Rankin J.M., 2001, MNRAS, 322, 438 
Franco L.M., Link B., Epstein R.I., 2000, ApJ, 543, 987 
Gil J., Melikidze G.I., Geppert U., 2003, A&A, 407, 315 
Goldreich P., Julian W.H., 1969, ApJ, 157, 869 
Goldstein H., 1980, Classical Mechanics. Adison- Wesley, Reading, 
MA 

Gould D.M., 1994, PhD thesis, Jodrell Bank Observatory 
Graham-Smith F., 2003, Reports of Progress in Physics, 66, 173 
Jones D.I., Andersson N., 2001, MNRAS, 324, 811 
Landau L.D., Lifshitz E.M., 1976, Mechanics. Pcrgamon Press, 
Oxford 

Link B., 2003, Phys. Rev. Lett., 91, 101101 

Link B., Cutler C, 2002, MNRAS, 336, 211 

Link B., Epstein R.I., 2001, ApJ, 556, 392 

Link B., Franco L.M., Epstein R.I., 1998, ApJ, 508, 838 

Lyne A.G., Pritchard R.S., Smith F.G, 1988, MNRAS, 233, 667 

McCulloch P.M., Hamilton P.A., McConnell D., King E.A., 1990, 

Nature, 346, 822 
Parry J., Hyslop H.R., Stairs I.H., Lyne A.G and Kramer M., 
2005, "Profile Changes in PSR B1828-11", Neutron Stars at 
the Crossroads of Fundamental Physics, display presentation 
Rankin J.M., 1990, ApJ, 352, 247 
Rankin J.M., 1993, ApJ, 405, 285 
Rankin J.M., Wright G.A.E., 20 03, AfcAR, 12, 43 
Ruderman M.A., 2001, preprint ( astro-ph/0109353 I 
Ruderman M.A., Sutherland P.G., 1975, ApJ, 196, 51 
Sedrakian A.D., 2005, Phys. Rev. D, 71, 083003 
Sedrakian A.D., Sedrakian D.M., 1995, ApJ, 447, 305 
Sedrakian A.D., Wasserman I., Cordes J.M., 1999, ApJ, 524, 341 
Shaham J., 1977, ApJ, 214, 251 
Shaham J., 1986, ApJ, 310, 708 

Stairs I.H., Lyne A.G., Shemar S.L., 2000, Nature, 406, 484 
Stairs I.H., Athanasiadis D., Kramer M., Lyne A.G, 2003, ASP 

Conf. Ser. Vol. vr 
Tannanbaum H., Gursky H., Kellog E.M., Levinston R., Schreier 

E., Giacconi R., 1972, ApJ, 300, L63 
Wasserman I., 2003, MNRAS, 341, 1020 
Weisberg J.M., Taylor J.H., 2002, ApJ, 576, 942 
Weisberg J.M., Romani R.W., Taylor J.H., 1989, ApJ, 347, 1030 
Wright G.A.E., 2003, MNRAS, 344, 1041 



ACKNOWLEDGMENTS 

We thank Ingrid Stairs for providing us with the data and 
for valuable discussion. This research is supported in part by 
NSF AST-0307273 (Cornell University), NSF AST-0098728 
(Montana State University) and IGPP 1222R from LANL. 



REFERENCES 

Alpar M.A., Langer S.A., Sauls J. A., 1984, ApJ, 282, 533 
Alpar M.A., Sauls J.A., 1988, ApJ, 327, 723 
Bisnovatyi-Kogan G.S., Mersov G.A., Sheffer E.K., 1990, SvA, 
34, 44 

Bisnovatyi-Kogan G.S., Kahabka P., 1993, A&A, 267, L43 
Blaskiewicz M., 1992, PhD thesis, Cornell University 
Bondi FL, Gold T., 1955, MNRAS, 115, 41 
Cordes J.M., 1993, ASP Conf. Ser. Vol. 36 
Cutler C, 2002, Phys. Rev. D, 66, 084025 
Cutler C, Ushomirsky G, Link B., 2003, ApJ, 588, 975 
D'Alessandro F., McCulloch P.M., 1997, MNRAS, 292, 879 
Davis L., Goldstein M., 1970, ApJ, 159, L81 



Precession of the Isolated Neutron Star PSR B1828-11 13 



APPENDIX A: PERIOD RESIDUALS FOR A TRIAXIAL RIGID STAR IN THE ABSENCE OF 
EXTERNAL TORQUES: THE GEOMETRIC CONTRIBUTION 

Euler's equation for a freely precessing rigid body is, 

^ + fixL = (Al) 

and can be solved analytically in terms of Jacobian elliptic functions (see Landau & Lifshitz, 1976). Using the principal axes 
(Ii ^ I2 $S ^3) as the basis for the body (rotating) frame, we can express the components of the angular momentum unit 
vector L as, 



Li = — Ai cn(r, k 2 ) , Ai 



'h(2EI a 


-I?) 


L 2 (h - 


-h) 


h{2Eh 


-L 2 ) 


L 2 (h - 


-h) 



Li = - A, si, ( r. k") , A 2 = yJ h L^*_j^ - = A i Vl + e 2 I ^ ' 

> h{L 2 -2Eh 
LHh-h) 

Here, the argument of the elliptic functions is 



L 3 = A 3 dn(r, k 2 ) , A 3 = x / 13 ^ { ~_j^ = yr^A? 



, (I 3 - h)(L 2 - 2Eh) eLA 3 
r = ft* where u p = ^ — = ] - ? _ (A3) 

and, the parameter of the elliptic functions is, 

2 _ (I 2 -h)(2Eh-L 2 ) _ e^A? _ 2 2 
k ~ {h~h){L 2 -2Eh) - A§ - eX (A4) 

where, we make use of the following auxiliary definitions, 

e={h-h)/h and e 2 = |||^] (AS) 

The minus signs that we have explicitly included in our definitions are due to our choice of the initial orientation of axes. 

Note that uj p is not the precession frequency, since the elliptic functions do not have a period of 2%, or more precisely, u) v 
is not the time derivative of the angular displacement. Instead, the precession frequency is given through, 

n P = ^ = ^ (A 6) 

where 2n is the period of the elliptic functions and can be calculated using the Legendre elliptic integral of the first kind, 
r = F(4>, k) where sin (f> = snr, 

tt/2 = F(tv/2, k) (A7) 

The values of the parameter k 2 are unrestricted, though different regimes require careful handling, k 2 < 1 corresponds 
to precession around the body z axis; k 2 — 1 corresponds to the unstable trajectories of the angular momentum, which decay 
exponentially towards the intermediate axis, y; and k 2 > 1 is precession around the x axis (see Binet ellipsoid). For now we 
will confine ourselves to the first case, and the other two will be left for a later section. 

We are interested in an isolated neutron star, and we want to determine the time of arrival (TOA) of pulses (produced 
along the magnetic axis) for an inertial observer. Let the inertial z axis be along the angular momentum vector, which 
remains constant; and let the inertial x axis be defined by the orientation of the observer, whom we choose to locate in the 
first quadrant of the inertial xz plane. Let a unit vector b denote the orientation of the magnetic axis, and bi be the rotating 
frame components. Then, whenever the inertial y component (which we choose to denote by by) vanishes, while the inertial x 
component (b x ) is positive, we get a pulse. The two frames are related through a rotation matrix constructed from the Euler 
angles 9,ip and <j> (see Goldstein 1980), whence the two conditions can be expressed as, 

by — 61 (cos V> sin <j> + cos 6 cos (psm ip) + 6 2 (— sin tp sin cj> + cos 6 cos ^>cos ip) — 63 sin 6 cos (f> — (A8) 

and, 

b x = bi(costj) cos 4> — cos 9 sin 4>simj)) — fe 2 (sin ip cos <f) + cos 9 sin 0cos tp) + 63 sin 9 sin <j> > (A9) 
Using the solution for the angular momentum (Eq. <A2I ). the Euler angles are given through, 

cos (9 = A3 dnr , sin# = Ai \J \ + e 2 sn 2 r , cos ip = — Vl + e snr ^ s j n ^ = cnr (A10) 

Vl + e 2 sn 2 r y/1 + e 2 sn 2 r 

dcp _ L_ . 
dt I 3 
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Figure Al. Orientation of the magnetic axis in the body frame. We assume that the pulses are also emitted along the same axis. 



The last equation can be written as, 

,/.-, , L <Jl + e 2 f T dr , 

m = ^o + T t+— — / 2 2 ah 

h A 3 J 1 + e 2 sn 2 r 

Eq. IA8I also implies, 

tan (p = — 



JV = fe 3 Ai (l + e 2 sn 2 r) + 6 2 A 3 snr dnrVl + e 2 + 6iA 3 cnrdnr (A12) 



D = 62 enr — 61 snry/ 1 + e 2 

Pulses are seen when Eqs. I)A11(I and l|A120 are both satisfied. In the absence of precession (Ai = and k 2 = 0) the solution 
of Eq. HA 1211 for <f> is simply given through, 



4> = 27m + 77 + tan ^VTTe 2 " tanr) (A13) 

where r\ = 7r/2 — </3 (Fig. lAlt : and we can further restrict it to lie anywhere between and 271". Note that we have implicitly 
included the second requirement, Eq. (IA9L by skipping every other possible solution for <j>. This is a non-trivial assumption, 
and would break down if the magnetic axis b happens to lie between tt and L. However, for a pulsar these two vectors are 
very nearly aligned since e is extremely small. Therefore, we do not need to worry about such a case. 
We express the general solution for <f> as, 

4> = 2nn + r) + C (A14) 
where £ is confined to lie within a period of tangent (i.e. n). Using tan 77 = &1/62 one can show that, 

tanC^-ff 1 (M5) 
s TVbi + Db 2 

We now have two equations for <j> (Eqs. lAlU and lAm ) which we can combine to get the times of arrival of pulses, 



L yl + e 2 1 dr , . , „. 

-tn^ + Cn-Co—^ TT72 - ?7 (A16) 



Co (for t = 0) appears as a consequence of the fact that <f> = r) + We have, 

b2baAi - 6162(1 - A :i ) 
6163A1 +61 + 6 2 A 3 

Note that in the absence of precession (i.e. when Ai = 0) the times of arrival reduce to the form, 

tn = ^ (A18) 

which is the solution for pure rotation. 

The period (between two pulses) is given as, 

P n = t„ - t n -i (A19) 
whence, 

L Pn -2n = ±AP n = C n -C«-i-^^ i,t 2 ( A2 °) 

h h A 3 y 1 + e 2 sn 2 r 
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If the precession period is much longer than the pulse period (as is the case for a neutron star) we can approximate the 
differences by derivatives, 



L a p _ dC - n Vl + eVAa dr„ 



dC, n 7T+^/A 3 

dr n 1 + e 2 sn 2 r„ 



73 dn 1 + e 2 sn 2 r n dn 

We will find it convenient to define the expression inside the parentheses as a new function, 



£ (A21) 



aY„ I + e 2 sn 2 T n 
The derivative d^/dr is given through (from Eq. IA15I 1 

rfC _ DN' - JW 
oY ~~ N 2 + D 2 

where, from Eq. 1A12I . 

dN 



u = ^ - --cr (A2 2) 



(A23) 



2 2 2 2 / 2 2 2 

= 2&3Aie snr cnr dm - + 02A3 cnr( dn r — A; sn r)vl + e 2 — &1A3 snr( dn r + cn r) 
ar 

dD 



hi snr dnr — 61 cnr dnr \f 1 + e 2 
or 

To evaluate dr n /dn we will make use of the time of arrival equation, Eq. (|A16|l . Taking the derivative of both sides with 
respect to n, we get, 

L dt n L dr„ dr n dr n 2tyzu p /aca\ 

— — = — = 2tt + /„— so that — = (A24) 

13 an IzLOp an an an 1 — w v j n 

where we have defined a new dimensionless quantity, 

--*?-7rh (A25) 

The pulse period will be given through, 

P n = t„ - tn-i ~ ^fc = — ^ = , P * (A26) 
an uj p an 1 — zu p j n 

whence the period residuals can be found to be, 

— = — - 1 = w,/,, (A27) 

P* P* l-Wp/n W ^ ^ 

where P* = 27r73 /L is the rotation period of the star, and the last approximation results from our anticipation that e will be 
sufficiently small. Indeed, from the above definitions we get, for small k 2 , 

" 3.2xlO- 8 ^i SeC j (A28) 



VTTe 2 P P ' P P (yrs) 

The coefficient of the function f n in Eq. I|A27(I is, 

B = P*vj p where vj p = — ^ = - — = — — (A29) 

L LP P irP p 

in other words, 

B=^ = ^ (A30) 
For PSR B1828-11 the rotation period is 405.04 ms, and the precession period is about 511 days, so that B ~ 3.8 ns. 



Al The Axisymmetric Body 

For an axisymmetric body we have e 2 = k 2 — 0. We can also set 62 = which is equivalent to introducing some initial phase 
shift in the definition of r. We thus get, after some rearrangement, 



fn 



d( 1 Ai &1&3A3 cost + ofAi sin r + &3A1 

dr A 3 A3 (63A1 + 61 A3 cos r) 2 + b\ sin 2 r 



Let's now assume that the angle 9 between the symmetry axis and the angular momentum is small, i.e. Ai 
A 3 ~ # 2 /2, and working to second order compute the period residuals, 



fn 



1 



+ 



cos 2r 



(A31) 
d and 

(A32) 
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Also let £13 = cos x- Then, 



1 



2 



APn/P* W U7 p f n « — #e cot x cos T — h# e — + cot x cos 2t (A33) 



(Here P* = 2-71-/3 /L and S7 P = IsLUp/L = eA3.) Note that a harmonic arises from geometrical effects. 

A2 Precession Around the Principal Axis Corresponding to the Smallest Moment of Inertia (k 2 > 1) 

We will now look at this case in more detail. The solution given through Eq. <A2> is still valid. However, it will be mathe- 
matically and computationally convenient to carry out a transformation of the Jacobian elliptic functions with k 2 > 1 into 
functions with parameter 1/fc < 1. (That, in the limit Ai = 1 and A3 = 0, we get uj p — whence r = and k = 00, provides 
further motivation.) Then, we can write the general solution as, 

L\ = — Ai dnf 



L 2 = - A 3 \/l + e 2 snf ( A34) 

L3 = A 3 cnf 

where, f = rfc, e 2 = 1/e 2 and the parameter of the elliptic functions is now k 2 = l/fc 2 . We can also define a new frequency 
from Eq. 1A3I . 



(I 2 - h){2Eh ~ L 2 ) 

uj p = ujpk = 4/ — — (A35) 

V I1I213 

Through an appropriate redefinition of axes, the solution can be expressed in a form identical to the k 2 < 1 case, except that 
now ii > I2 > J3. Define a new right-handed coordinate basis for the rotating frame, 

ei = — e.3 

e 2 = -e 2 (A36) 
§3 = — ei 

Let Ai = A3 and A3 = Ai , then the components of the angular momentum can be expressed as, 
L\ — —L3 = — Ai cnf 



L 2 = -L 2 = AiVl +e 2 snf (A37) 
L3 = —Li — A3 dnf 

The precession is now clockwise, as can also be verified from Euler's equation. Ai have exactly the same form as before, in 
terms of the new moments of inertia, 



A 3 = 



'h(L 2 - 


2Eh) 


L 2 (h 


-Is) 


l h{L 2 ^ 


2Eh) 


L 2 (i 2 


-h) 


l h{2Eh 


-L 2 ) 



Ai 



A * = \^77- f TT^ (A38) 



L 2 (h-h) 



So do tip, e and k , as can be verified from the equations above. We have thus transformed the problem from a "k > 1 case 
for an 73 > ii body" into a "k 2 < 1 case for an 7i > 73 body" , which should not be surprising. 

The equations for the Euler angles (Eqs. iAIOi ~) remain of the same form, with the exception of cosip. This is effectively 
a sign change, r — > — f , in the argument, 

C(f) = C(-f) whence ^ = -^-^ (A39) 
One must be careful with Eq. (IA1H as well, where there is also a sign change due to the fact that now 7i > 73, 



Hh^hl = _VTTP = __A^ (A4Q) 

ujphh A 3 A!A 3 

These two effects add up to modify the function /„ defined through Eq. 1A22L 

fn(r)=-fn(-f) (A41) 
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APPENDIX B: THE WOBBLE AND BEAM SWING ANGLES 

We will define the wobble (9) and beam swing {•&) angles according to (see Fig. , 

cos# = L z = A 3 dnr (Bl) 
cos ■& = Lb = — 61 Aj cnr — b 2 A 2 snr + &3A3 dnr (B2) 

By definition, the wobble angle is equivalent to Euler's angle 6 and is constant for an axisymmetric star. The beam swing 
angle is related to the angle between the beam and the observer. It could exceed 90°, but since the pulse will have a limited 
angular size, there is a restriction on how much it can vary throughout a precession period. Otherwise, the beam will leave 
the observer's line of sight. Therefore, the span of the beam swing angle serves as a constraint. The angle can be further 
restricted by imposing the conditions for an interpulse. 

We define the widest span as Ai? = # ma x — $msn- Then the constraint is that this be smaller than some value A'&max, which 
is estimated based on information about the pulse width and shape. Note that the beam angle depends on four parameters: 
the two angles determining the orientation of the pulse, and any two of k 2 , e 2 and A = A1/A3. 



APPENDIX C: PERIOD RESIDUALS FOR THE SPINDOWN TORQUE 

When external torques are present Euler's equation becomes, 

^ + fixL = N (CI) 
Taking the dot product with the angular momentum, we get the equation governing the evolution of its magnitude, 

f=L-N (02) 
If we now substitute L = LL in Euler's equation, we get, after some rearrangement, 

L^t +L 2 (I- 1 L) x L = N- (L- N)L (03) 

which governs the evolution of the orientation of the angular momentum. These two are the basic equations that need to be 
solved. Of course, only three (of the total of four components) are independent equations. 

In the classical rotating magnetic-dipole model of pulsars, the angular momentum is lost to radiation. The electromagnetic 
torque for a spherical rigid star in vacuum is (Davis and Goldstein, 1970), 

N„ ac = -^^bx(nxb) = -^^[f2-(n.b)bl (04) 

6c Li 3c d L J 

Note that the torque vanishes when ft and b are aligned. However, the pulsar is not in a perfect vacuum, and is surrounded 
by a magnetosphere. Therefore, there should be loss of angular momentum even when these two vectors are aligned. We will 
therefore adopt a general spindown torque of the form, 

N sd = —N [n - o(n ■ b)b] (C5) 

where a is a dimensionless parameter that measures the relative importance of the two components. The amplitude of the 
spindown torque can be estimated from observed values of the spindown time, and is small. We will be interested in a particular 
example (PSR B1828-11) where the spindown time is, 

tsd ~ jjr ~ 10 5 yrs (06) 
Compare this with the observed precession period for the same pulsar, 

P, = gL,* lyr (07) 
The ratio of the two gives, 

(C8) 

The second term in Eq. ()C3fl has a magnitude of eL 2 /I, therefore the RHS of that equation is quite negligible for the case of 
interest (as will be discussed below). 

The above form of the torque is true for spherical stars. This is nevertheless a good approximation, given how small N 
and e are. In fact, we will neglect all combinations of N with e. This is equivalent to taking f2 ~ L within all torque terms, 
since the angle between these two vectors is of the order of e. We therefore have, 

N = -N„ [t - a(t ■ b)b] (09) 
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^ = L ■ N = -N [l - a(L • b) 2 ] (CIO) 

L^+i 2 (r'L)xL = N-(L-N)L = JV a(L ■ b) lb - (L • b)L| (Oil) 
at 1 J 

Loss of energy (or angular momentum, given through the second equation) now clearly requires that a ^ 1. Finally, we will 
also neglect any time dependence within N itself. 

The third equation demands careful thought. In component form, 

(bi — L\ cos ■& \ 
62 -£2 cos (C12) 
fa — L s cos # / 

where Li are the components of L, s — (I2 — — Ii), cosi? = L ■ b, and we have already neglected second order terms 

in e. As long as the star is sufficiently non-spherical and the angular momentum is sufficiently misaligned with the body z 
axis, we can neglect the RHS, as it causes changes in the orientation of the angular momentum smaller (by many orders of 
magnitude) than the second term. In other words, e and Ai are small but not zero. (Keep in mind that there is no precession 
if either one is zero.) Also, s cannot be too close to unity, i.e. e 2 cannot be exceptionally large. 

The same cannot be done in the RHS of the equation for the magnitude of the angular momentum, as it is the only term 
we have. Incidentally, setting N — would take us back to the torque-free precession case. 

In order to write the equations in a dimensionless form, let's divide all sides by a frequency lu p , defined in accordance 
with Eq. llAll . 

Mt) = (ci3) 

but where the magnitude of the angular momentum is no longer constant. Also define, 

dr = Lu p (t)dt (C14) 

which, for a constant lu p , reduces to the familiar form of the torque-free case. Now, the differential equations become, after 
some rearrangement, 

^ = L-N/o;,, = -(JVcMO [l -a(L- b) 2 ] (C15) 

^ + (%)(r 1 L)xL = (C16) 
ar 

Since L/uj p is time-independent, the second equation has exactly the same solution as before, except that r is now different, 
and given through a differential equation on its own. In other words, Li's remain of the same form. Thus, we only need to 
solve Eqs. OH and (IHTst . 

Define a new dimensionless constant, 



vj p = hu v /L = eA 3 /Vl + e 2 < e (C17) 
and let's write, 

L = L [1-£(t)] 

N = -A D n(r) (C18) 

t = [T + S(t)]/lU P o 

where u po — zu p L /Is. The differential equations now become, 
dl ( hN \L-n dS £ 

Tt = Y^i? — and Tr^Y-i (C19) 



It's worth noting that we make no assumptions in these substitutions 
If we finally define one more dimensionless constant, 

zu p L% 



r sd = (C20) 



and let I = Y s di and 8 — T s dS, then the two equations can be written as, 

# = ^-_ and f = — U (C21) 
For the pulsar that we discuss here, we have, 

* hN P p 10 _ 5 



tL% t sd 



(C22) 
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This means that one may safely ignore the denominators of the two equations, thus further simplifying the results, 

4^=L-n and ^- =i (C23) 

dr dr 

Note that L n = 1 — a(L-b) 2 ^ 0, thus assuring that L is monotonically decreasing (as required by loss of angular momentum). 



CI Time of Arrival Residuals 

Since Li remain of the same form, with the only difference being that r is now determined through a differential equation 
(Eq. IIC 1411 ), the Euler angles remain the same (Eqs. IA10I 1. However, due to the time dependence of L, it is more convenient 
to express as a function of r, and we need to replace Eq. JA1 111 by, 



b + 



dr 



l 



b + + — 



dr 



1 + e 2 sn 2 r 



(C24) 



Thus, all we have to do is to replace Lt n /I^ by r n /zu p on the LHS of Eq. JA16I I. The period (given through Eq. l|A26)l 1 
remains the same as well, as does the calculation of dr n /dn. In fact, we run into trouble only with the period residuals, since 
the magnitude of the angular momentum is now changing. Define, 

2tt/ 3 , \ n T~k 27T/ 3 



APres — Pn 



L 



and 



AP n = P n 



(C25) 



AP n is formally the same as the torque-free case. However, observations give us only information about AP res . In practice, 
one first determines the period (P*) at some epoch (t a ), and then finds the period derivative (P*) which is the secular term 
attributed to spindown, and subtracts both contributions, so that the residuals are then given through, 



AP res = P(t) - P* - P* (t - U) 
Consider the difference between the two definitions above, 



AP„ 



AP n 



2./a (I 



1 

T 



2-Klst a d dS 



L dr 

where we make use of Eqs. l|cT9t and gS). We thus get, 



AP r , 



1 - VJ p f n 1 - f sd 



(C26) 



(C27) 



(C28) 



where P* = 2n/Q+ = 2nIs/L . The first term is the geometric effect (AP ge /P* w vj v f n ) and the second term is the spindown 
term (AP a d/P+ w T a d£). There are still secular terms present in the spindown term that need to be subtracted. This will be 
taken care of below. The relative amplitude of these two terms cannot be simply determined, and it is possible that either 
one is dominant, or that they are comparable. 

We now turn our attention to the calculation of £. From Eq. 1C23I we have, 



1 = 



L • n dr 



where, 



L-n = l-a(L-b) 2 



and, 



L ■ b = — 61A1 enr — 62A2 snr + 63 A3 dnr = cos x 



(C29) 



(C30) 



Carrying out the integrals of the Jacobian elliptic functions, and substituting the values of k and A2, we get, after some 
rearrangement, 



(L- 



b 2 2 (l + e 2 



1.2 ; 2 



T + 



, 2 



+ e) + b 2 3 e 2 



E(am t, k) 



(C31) 



+ 



26i6 2 A| v / T+i 5 



(1- dnr)-2AiA 3 6263 VI + e 2 (l — enr) + b\bs si 



where we have introduced the complementary parameter k\ — 1 — k 2 . P(amr, k) is the Legendre elliptic integral of the second 
kind, and amr is the Jacobi amplitude, amr = sin -1 snr. Note that we have implicitly assumed that at the zero of time, 
the angular momentum is in the xz plane of the body frame. This is a non-trivial assumption, and in general one does not 
have the freedom of randomly setting the initial orientation of the angular momentum. Therefore, in general we would have 
L = L(r — t ), where t is the phase offset and is an additional parameter, and the overall result would be to replace 
above by i[j — r ) — £(— t ), where £(— r D ) is just a constant. For simplicity, we will continue to assume r = in the rest of 
our derivations, but the general case should be kept in mind. 

The oscillatory part of I, after secular terms have been removed, is given through, 



A£: 



(L • n)r 



(C32) 
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The average is carried over a precession period, 2tt = 4F(ir/2, k), 

(L- n) = 1 -a((L-b) 2 } 
where, using Eq. fCTH we get, 



<(L-bf> = - 



r A 2 
(L-b) 2 dr = ^ 



b 2 2 (l + e 2 ) - b 2 k 2 



b{-bl{l + e 2 ) + ble 2 



E(*/2,k) 

F(7T/2,fc) 



and we have made use of the relations am(27r) = 2tv and E(2ir, k) = 4E(-k/2, k). We thus get, 
Al = a{(t ■ b) 2 )r - a [ (L • b) 2 dr 



which can now be used in Eq. I|C280 to calculate the time of arrival residuals, 
^ « f sd Ai 

We will find it convenient to express this equation in the following form, 



Ai/a = ci(l — cnr) + c 2 snr + — dnr) + p- 



B(tt/2,*) 



t — E(&m t, k) 



where, 

ci = 2A2A36263 , c 2 = 2A1A3&1&3 , C3 
These coefficients are related to each other through 

C4 



F(7T/2,fe) 

2AiA 2 fei6 2 and c 4 = A? [61 - fe 2 (l + e) + bie z 



(C33) 
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It is also interesting to note that A£ has a non-zero average over a precession period. The residuals may have non-zero 
average depending on when and how the period and its derivatives are calculated. This becomes particularly important when 
calculating the time of arrival residuals, which can be obtained by integrating the period residuals, and if the period residuals 
have a constant term, then the time of arrival residuals will have a linear term. Therefore, in calculating the time of arrival 
residuals one will have to subtract any constant terms from the period residuals. 

C2 Amplitude of the Residuals 

Consider the period derivative, which is given through the secular terms in £, 

P* = rjP^T 3d u p (C39) 
where r\ = (L • n) = 1 — ac and c = (cos 2 d). The amplitude of the period residuals thus becomes, from Eq. (ICJ36I . 



aP+ 



aA a 

lu p (1 — ac a ) 1 — ac a 



where A a 



2io v t s 



(C40) 



and t s d — P+/2P+ is the spindown time. Recall that a measures the strength of the oscillating part of the spindown torque, 
and must be ^ 1. 

For PSR B1828-11 the period is 405.04 ms, the precession period is about 511 days, and the spindown time is 0.11 Myr, 
so that we get A a ~ 409.95 ns. 



C3 The Axisymmetric Body 

For an axisymmetric star e 2 = k 2 — 0, but A = k/e 7^ 0. Due to the symmetry we can set 6 2 = by shifting the zero of time 
through some phase r . (Note that the same cannot be done in a triaxial body, where we chose to fix the axes according to 
the principal moments of inertia.) In this case Eq. l|C37fl reduces to the form, 



M/a 



A 



1 + A 2 



sin 2\ sin(r — r ) — — sin 2 \ sm 2(t — r Q ) 
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Note that the non-linearity of the dipole contribution of the torque naturally brings in a harmonic. The period residuals are 
then, 



AP sl 

Po 



T s dAl where 



N 
ujpL 



Here r c is the characteristic time, 
_ 3c 3 J 3 



(C42) 
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It turns out that, for the axisymmetric case, the geometric term is quite negligible compared to the spindown term, for the 
range of physical parameters of interest (Jones & Andersson, 2001; Link & Epstein, 2001). 

To convert our result for period residuals (AP/P ) into residuals of the derivative of the angular velocity (Afl/fl ) given 
by Link & Epstein, we make use of, 

. + dAP dAP . Ail AP 

AP=——=Up— — and — — = — (C44) 

dt P dr n o P K ' 

which indeed gives the correct results, together with the initial phase difference of 7r between the two definitions, 
Af2 a\ 



o 2r c (l + A 2 ) 



A 2 

- sin 2xcos(r — r ) + — sin x cos 2(t — r Q ) 
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There is one difference between the two derivations and that is the presence of the coefficient a which measures the 
strength of the spindown torque. With the addition of this new element, the number of unknowns increases to three (a, A and 
X), while a fit to data will yield only two coefficients (di and t does not contain any further information). This implies 
that there is a certain level of freedom in the choice of the physical coefficients. Let's denote the fitting function by /, 

f = ai sin(r — r ) — £12 sin 2(r — To) (C46) 

Then, the relations between the coefficients of this function and the physical parameters that we actually seek would be, 

aAsin2y , aA 2 sin 2 y ,^,„n 

ai = TTA^ and a 2 = WT ^ (C47) 

Define A = tan 9, and the ratio of the two coefficients gives, 

tanxtan(9 = — (C48) 

It is also possible to express x an d as functions of a. However, as it turns out, the range of the physical parameters is 
severely restricted by the beam swing angle constraint, which in the axisymmetric case is given through, 

Atf = 2 min ( X , 0) < Ai? mM (C49) 

This forces one of the two angles to be small (which will have tani < 0.09 even if we let Afi max = 10°); while Eq. iCl8l 
ensures that the other remains very close to 90°. (The ratio of the coefficients is found to be 02/01 ~ 0.4 for the data used 
by Link & Epstein (2001). This yields the condition tan2 > 36, i.e. the second angle has to be larger than 88°, in accordance 
with previous findings.) 



APPENDIX D: STATISTICAL INFERENCE 

Denote the set of parameters by x. Then the pdf for the parameters can be calculated as, by Bayes's theorem, 
P{x\M)P(D\x,M) 



P(x\D,M) = 



P{D\M) 



(Dl) 



where D stands for data, M stands for the model and also takes into account any other information that is available on 
the problem apart from data (in this case the beam swing angle, which we impose as a restriction on the parameter space). 
P(x\M) is the prior probability for the parameters; P(D\x, M) is the likelihood; and P(D\M) is effectively a normalization 
constant. 

To find the pdf for a certain parameter, or a subset of parameters, we integrate Eq. IDlll over the remaining parameters. 
For that we need to know the likelihood. For each data point yi at time ti, we have a theoretical prediction /, = f(ti\x, M). 
For well-known uncertainties with Gaussian distribution, we would then have, 



P(D\x, M) = JJ^V^r 1 exp 



(Vi - 



2a? 



(D2) 



If we assume that the error bars are not well-determined and rescale them through some number F, the above equation 
becomes, 



P(D\x, M) = JJtFcriV^)- 1 exp 



2F 2 c 



(D3) 



and we regard F as an additional parameter. We have to introduce a prior for F. Since we do not want it to depend much on 
the endpoints, we take it to be flat over d In F, i.e. proportional to dF/F, and integrate over all values of F. Define, 



E 



(Vi - fi) 2 
2a? 
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whence we get, for d data points, 



f 

Jo 



dF P(D\x,M) (tt 1 \ f 00 e-^ p2 dF ra^-i-d/2 m _. 



where we have dropped anything that does not depend on the remaining parameters x, including integrals that give constants, 
products of the original sigmas, and factors of \/2%. Thus, our final result is, 

P(x\D,M) oc P(x\M) [x 2 o(x)]~ d/2 (D6) 

where the first term is the prior probability, and the constant of proportionality can be computed from the condition that the 
final pdf is normalized to one. 



